clear; close all;
saveArg = 0;


%% Plotting

medSize = 70;
w = 2;
MadMedTheshold = 0.6;

labelSzAx = 14;
labelSzXSec = 14;
labelPlot = 18;
labelLeg = 12;
labelDeviations = 12;

xPos1 = 0.13;   xPos2 = 0.53;
yPos1 = 0.72;   yPos2 = 0.42;   yPos3 = 0.12;
wPos1 = 0.21;   wPos2 = 0.35;   hPos = 0.23;
xleg = 0.88;    yleg = 0.58;    wleg = 0.09;    hleg = 0.13;

figure(1000); clf;
set(gcf,'units','normalized','OuterPosition',[0.2 0.05 0.5 0.95]);





% Circulation
data = importdata('Exp_21OCTSolid\Data_03-15000-15100\B00001.dat');
data = data.data;
x = data(:,1);
y = data(:,2);
vx = data(:,3);
vy = data(:,4);
noInd = find(vx==0&vy==0);
x(noInd) = [];
y(noInd) = [];
vx(noInd) = [];
vy(noInd) = [];
vector = sqrt(vx.^2+vy.^2);
theta = acos(vx./vector);
threshold = max([10^-5,median(vector)-mad(vector)]);

xlin = linspace(min(x),max(x),medSize);
ylin = linspace(min(y),max(y),medSize);
uMed = nan(medSize,medSize);
vMed = uMed;
uMad = uMed;
vMad = uMed;
for i = 2:medSize-1
for j = 2:medSize-1
    indCheck = x>xlin(i-1) & x<xlin(i+1) & ...
        y>ylin(j-1) & y<ylin(j+1) & vector>threshold;
    uMed(j,i) = median(vx(indCheck));
    vMed(j,i) = median(vy(indCheck));
    uMad(j,i) = mad(vx(indCheck));
    vMad(j,i) = mad(vy(indCheck));
end
end
uMed(abs(uMad)./abs(uMed)>MadMedTheshold) = NaN;
vMed(abs(vMad)./abs(vMed)>MadMedTheshold) = NaN;
reducedVector = sqrt(uMed.^2+vMed.^2);

thetaMed = acos(uMed./reducedVector);
[xMed,yMed] = meshgrid(xlin,ylin);

thetaMed = reshape(thetaMed,numel(thetaMed),1);
xMed = reshape(xMed,numel(xMed),1);
yMed = reshape(yMed,numel(yMed),1);


indA = xMed>16-w & xMed<16+w;
xa = xMed(indA);   ya = yMed(indA);   za = thetaMed(indA);
indB = xMed>-9-w & xMed<-9+w;
xb = xMed(indB);   yb = yMed(indB);   zb = thetaMed(indB);
indC = yMed>-233-w & yMed<-233+w;
xc = xMed(indC);   yc = yMed(indC);   zc = thetaMed(indC);
indD = yMed>-139-w & yMed<-139+w;
xd = xMed(indD);   yd = yMed(indD);   zd = thetaMed(indD);


quivSkip = 7;
indX = unique(x);
indX = indX(1:quivSkip:end);
indY = unique(y);
indY = indY(1:quivSkip:end);
indFlagX = zeros(size(x));
indFlagY = indFlagX;
for i = 1:length(indX)
    indFlagX(x==indX(i)) = 1;
    indFlagY(y==indY(i)) = 1;
end
indFlag = indFlagX==1 & indFlagY==1;

subplot(3,2,1);
scatter(xMed,yMed,30,thetaMed,'filled'); c=colorbar; hold on;
% quivSkip = 10;
quiver(x(indFlag),y(indFlag),vx(indFlag),vy(indFlag),6,'k');
plot([xa(1),xa(end)],[ya(1),ya(end)],'c-','linewidth',3);
plot([xb(1),xb(end)],[yb(1),yb(end)],'c--','linewidth',3);
plot([xc(1),xc(end)],[yc(1),yc(end)],'k-','linewidth',3);
plot([xd(1),xd(end)],[yd(1),yd(end)],'k--','linewidth',3);
set(gca,'Fontsize',labelSzAx);
xlim([-100,150]);
ylim([-300,-50]);
% ylabel('y (mm)');
set(gca,'XTick',[]);
caxis([0,pi]); set(c,'Ticks',[0,pi],'TickLabels',[0,180]);
tx = text(14,-285,'A'); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(14,-60,'A'''); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(-17,-290,'B'); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(-17,-60,'B'''); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(-90,-230,'C'); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(125,-230,'C'''); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(-90,-135,'D'); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(125,-135,'D'''); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(-90,-70,'a'); set(tx,'Fontsize',labelPlot);

[yasort,indsort] = sort(ya);
zasort = za(indsort);
yasort = (yasort-min(yasort))/range(yasort);
yasort = yasort(isnan(zasort)==0);
zasort = zasort(isnan(zasort)==0);
[ybsort,indsort] = sort(yb);
zbsort = zb(indsort);
ybsort = (ybsort-min(ybsort))/range(ybsort);
ybsort = ybsort(isnan(zbsort)==0);
zbsort = zbsort(isnan(zbsort)==0);
[xcsort,indsort] = sort(xc);
zcsort = zc(indsort);
xcsort = (xcsort-min(xcsort))/range(xcsort);
xcsort = xcsort(isnan(zcsort)==0);
zcsort = zcsort(isnan(zcsort)==0);
[xdsort,indsort] = sort(xd);
zdsort = zd(indsort);
xdsort = (xdsort-min(xdsort))/range(xdsort);
xdsort = xdsort(isnan(zdsort)==0);
zdsort = zdsort(isnan(zdsort)==0);
std1a = mad(zasort)*180/pi;
std1b = mad(zbsort)*180/pi;
std1c = mad(zcsort)*180/pi;
std1d = mad(zdsort)*180/pi;


subplot(3,2,2);
plot(yasort,zasort,'c-'); hold on;
plot(ybsort,zbsort,'c--');
plot(xcsort,zcsort,'k-');
plot(xdsort,zdsort,'k--');
set(gca,'Fontsize',labelSzAx);
ylim([0,pi*1.7]);
set(gca,'YTick',[0,pi],'YTickLabel',[0,180]);
% ylabe5l('\theta (deg)');
set(gca,'XTick',[]);
chld = get(gca,'Children');
set(chld,'Linewidth',2);
tx = text(0.18,pi*1.33,{...
    sprintf('\\sigma_{A-A''} = %0.0f^{\\circ}   \\sigma_{B-B''} = %0.0f^{\\circ}',std1a,std1b);...
    sprintf('\\sigma_{C-C''} = %0.0f^{\\circ}   \\sigma_{D-D''} = %0.0f^{\\circ}',std1c,std1d)});
    set(tx,'Fontsize',labelDeviations);
tx = text(0.05,pi*1.5,'d'); set(tx,'Fontsize',labelPlot);





% Buoyant
data = importdata('Exp_Prev_C1\Data_24-10200-10300\B00001.dat');
data = data.data;
x = data(:,1);
y = data(:,2);
vx = data(:,3);
vy = data(:,4);
noInd = find(vx==0&vy==0);
x(noInd) = [];
y(noInd) = [];
vx(noInd) = [];
vy(noInd) = [];
vector = sqrt(vx.^2+vy.^2);
theta = acos(vx./vector);
threshold = max([10^-5,median(vector)-mad(vector)]);

xlin = linspace(min(x),max(x),medSize);
ylin = linspace(min(y),max(y),medSize);
uMed = nan(medSize,medSize);
vMed = uMed;
uMad = uMed;
vMad = uMed;
for i = 2:medSize-1
for j = 2:medSize-1
    indCheck = x>xlin(i-1) & x<xlin(i+1) & ...
        y>ylin(j-1) & y<ylin(j+1) & vector>threshold;
    uMed(j,i) = median(vx(indCheck));
    vMed(j,i) = median(vy(indCheck));
    uMad(j,i) = mad(vx(indCheck));
    vMad(j,i) = mad(vy(indCheck));
end
end
uMed(abs(uMad)./abs(uMed)>MadMedTheshold) = NaN;
vMed(abs(vMad)./abs(vMed)>MadMedTheshold) = NaN;
reducedVector = sqrt(uMed.^2+vMed.^2);

thetaMed = acos(uMed./reducedVector);
[xMed,yMed] = meshgrid(xlin,ylin);

thetaMed = reshape(thetaMed,numel(thetaMed),1);
xMed = reshape(xMed,numel(xMed),1);
yMed = reshape(yMed,numel(yMed),1);

indA = xMed>18-w & xMed<18+w;
xa = xMed(indA);   ya = yMed(indA);   za = thetaMed(indA);
indB = xMed>-31-w & xMed<-31+w;
xb = xMed(indB);   yb = yMed(indB);   zb = thetaMed(indB);
indC = yMed>-70-w & yMed<-70+w;
xc = xMed(indC);   yc = yMed(indC);   zc = thetaMed(indC);
indD = yMed>30-w & yMed<30+w;
xd = xMed(indD);   yd = yMed(indD);   zd = thetaMed(indD);


quivSkip = 7;
indX = unique(x);
indX = indX(1:quivSkip:end);
indY = unique(y);
indY = indY(1:quivSkip:end);
indFlagX = zeros(size(x));
indFlagY = indFlagX;
for i = 1:length(indX)
    indFlagX(x==indX(i)) = 1;
end
for i = 1:length(indY)
    indFlagY(y==indY(i)) = 1;
end
indFlag = indFlagX==1 & indFlagY==1;

subplot(3,2,3);
scatter(xMed,yMed,30,thetaMed,'filled'); c=colorbar; hold on;
% quivSkip = 20;
quiver(x(indFlag),y(indFlag),vx(indFlag),vy(indFlag),5,'k');
plot([xa(1),xa(end)],[ya(1),ya(end)],'c-','linewidth',3);
plot([xb(1),xb(end)],[yb(1),yb(end)],'c--','linewidth',3);
plot([xc(1),xc(end)],[yc(1),yc(end)],'k-','linewidth',3);
plot([xd(1),xd(end)],[yd(1),yd(end)],'k--','linewidth',3);
set(gca,'Fontsize',labelSzAx);
ylim([-300,150]);
xlim([-100,150]);
ylabel('y (mm)');
set(gca,'XTick',[],'YTick',-300:200:100);
caxis([0,pi]); set(c,'Ticks',[0,pi],'TickLabels',[0,180]);
tx = text(13,-300,'A'); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(13,140,'A'''); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(-37,-300,'B'); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(-37,140,'B'''); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(-90,-80,'C'); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(110,-80,'C'''); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(-90,30,'D'); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(110,30,'D'''); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(-90,110,'b'); set(tx,'Fontsize',labelPlot);

[yasort,indsort] = sort(ya);
zasort = za(indsort);
yasort = (yasort-min(yasort))/range(yasort);
yasort = yasort(isnan(zasort)==0);
zasort = zasort(isnan(zasort)==0);
[ybsort,indsort] = sort(yb);
zbsort = zb(indsort);
ybsort = (ybsort-min(ybsort))/range(ybsort);
ybsort = ybsort(isnan(zbsort)==0);
zbsort = zbsort(isnan(zbsort)==0);
[xcsort,indsort] = sort(xc);
zcsort = zc(indsort);
xcsort = (xcsort-min(xcsort))/range(xcsort);
xcsort = xcsort(isnan(zcsort)==0);
zcsort = zcsort(isnan(zcsort)==0);
[xdsort,indsort] = sort(xd);
zdsort = zd(indsort);
xdsort = (xdsort-min(xdsort))/range(xdsort);
xdsort = xdsort(isnan(zdsort)==0);
zdsort = zdsort(isnan(zdsort)==0);
std2a = mad(zasort)*180/pi;
std2b = mad(zbsort)*180/pi;
std2c = mad(zcsort)*180/pi;
std2d = mad(zdsort)*180/pi;


subplot(3,2,4);
plot(yasort,zasort,'c-'); hold on;
plot(ybsort,zbsort,'c--');
plot(xcsort,zcsort,'k-');
plot(xdsort,zdsort,'k--');
set(gca,'Fontsize',labelSzAx);
ylim([0,pi*1.7]);
set(gca,'XTick',[]);
set(gca,'YTick',[0,pi],'YTickLabel',[0,180]);
ylabel('\theta (deg)');
chld = get(gca,'Children');
set(chld,'Linewidth',2);
tx = text(0.18,pi*1.33,{...
    sprintf('\\sigma_{A-A''} = %0.0f^{\\circ}    \\sigma_{B-B''} = %0.0f^{\\circ}',std2a,std2b);...
    sprintf('\\sigma_{C-C''} = %0.0f^{\\circ}    \\sigma_{D-D''} = %0.0f^{\\circ}',std2c,std2d)});
    set(tx,'Fontsize',labelDeviations);
tx = text(0.05,pi*1.5,'e'); set(tx,'Fontsize',labelPlot);







% Erupting
data = importdata('Exp_Prev_C1\Data_27-11400-11500\B00001.dat');
data = data.data;
x = data(:,1);
y = data(:,2);
vx = data(:,3);
vy = data(:,4);
noInd = find(vx==0&vy==0);
x(noInd) = [];
y(noInd) = [];
vx(noInd) = [];
vy(noInd) = [];
vector = sqrt(vx.^2+vy.^2);
theta = acos(vx./vector);
threshold = max([10^-5,median(vector)-mad(vector)]);

xlin = linspace(min(x),max(x),medSize);
ylin = linspace(min(y),max(y),medSize);
uMed = nan(medSize,medSize);
vMed = uMed;
uMad = uMed;
vMad = uMed;
for i = 2:medSize-1
for j = 2:medSize-1
    indCheck = x>xlin(i-1) & x<xlin(i+1) & ...
        y>ylin(j-1) & y<ylin(j+1) & vector>threshold;
    uMed(j,i) = median(vx(indCheck));
    vMed(j,i) = median(vy(indCheck));
    uMad(j,i) = mad(vx(indCheck));
    vMad(j,i) = mad(vy(indCheck));
end
end
uMed(abs(uMad)./abs(uMed)>MadMedTheshold) = NaN;
vMed(abs(vMad)./abs(vMed)>MadMedTheshold) = NaN;
reducedVector = sqrt(uMed.^2+vMed.^2);

thetaMed = acos(uMed./reducedVector);
[xMed,yMed] = meshgrid(xlin,ylin);

thetaMed = reshape(thetaMed,numel(thetaMed),1);
xMed = reshape(xMed,numel(xMed),1);
yMed = reshape(yMed,numel(yMed),1);

indA = xMed>18-w & xMed<18+w;
xa = xMed(indA);   ya = yMed(indA);   za = thetaMed(indA);
indB = xMed>-31-w & xMed<-31+w;
xb = xMed(indB);   yb = yMed(indB);   zb = thetaMed(indB);
indC = yMed>-70-w & yMed<-70+w;
xc = xMed(indC);   yc = yMed(indC);   zc = thetaMed(indC);
indD = yMed>30-w & yMed<30+w;
xd = xMed(indD);   yd = yMed(indD);   zd = thetaMed(indD);


quivSkip = 7;
indX = unique(x);
indX = indX(1:quivSkip:end);
indY = unique(y);
indY = indY(1:quivSkip:end);
indFlagX = zeros(size(x));
indFlagY = indFlagX;
for i = 1:length(indX)
    indFlagX(x==indX(i)) = 1;
end
for i = 1:length(indY)
    indFlagY(y==indY(i)) = 1;
end
indFlag = indFlagX==1 & indFlagY==1;

subplot(3,2,5);
scatter(xMed,yMed,30,thetaMed,'filled'); c=colorbar; hold on;
% quivSkip = 10;
quiver(x(indFlag),y(indFlag),vx(indFlag),vy(indFlag),5,'k');
plot([xa(1),xa(end)],[ya(1),ya(end)],'c-','linewidth',3);
plot([xb(1),xb(end)],[yb(1),yb(end)],'c--','linewidth',3);
plot([xc(1),xc(end)],[yc(1),yc(end)],'k-','linewidth',3);
plot([xd(1),xd(end)],[yd(1),yd(end)],'k--','linewidth',3);
set(gca,'Fontsize',labelSzAx);
ylim([-300,150]);
xlim([-100,150]);
set(gca,'YTick',-100:200:300);
xlabel('x (mm)');
% ylabel('y (mm)');
caxis([0,pi]); set(c,'Ticks',[0,pi],'TickLabels',[0,180]);
tx = text(13,-300,'A'); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(13,140,'A'''); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(-37,-300,'B'); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(-37,140,'B'''); set(tx,'Fontsize',labelSzXSec,'Color','c');
tx = text(-95,-80,'C'); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(105,-80,'C'''); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(-95,30,'D'); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(105,30,'D'''); set(tx,'Fontsize',labelSzXSec,'Color','k');
tx = text(-90,130,'c'); set(tx,'Fontsize',labelPlot);

[yasort,indsort] = sort(ya);
zasort = za(indsort);
yasort = (yasort-min(yasort))/range(yasort);
yasort = yasort(isnan(zasort)==0);
zasort = zasort(isnan(zasort)==0);
[ybsort,indsort] = sort(yb);
zbsort = zb(indsort);
ybsort = (ybsort-min(ybsort))/range(ybsort);
ybsort = ybsort(isnan(zbsort)==0);
zbsort = zbsort(isnan(zbsort)==0);
[xcsort,indsort] = sort(xc);
zcsort = zc(indsort);
xcsort = (xcsort-min(xcsort))/range(xcsort);
xcsort = xcsort(isnan(zcsort)==0);
zcsort = zcsort(isnan(zcsort)==0);
[xdsort,indsort] = sort(xd);
zdsort = zd(indsort);
xdsort = (xdsort-min(xdsort))/range(xdsort);
xdsort = xdsort(isnan(zdsort)==0);
zdsort = zdsort(isnan(zdsort)==0);
std3a = mad(zasort)*180/pi;
std3b = mad(zbsort)*180/pi;
std3c = mad(zcsort)*180/pi;
std3d = mad(zdsort)*180/pi;

subplot(3,2,6);
plot(yasort,zasort,'c-'); hold on;
plot(ybsort,zbsort,'c--');
plot(xcsort,zcsort,'k-');
plot(xdsort,zdsort,'k--');
set(gca,'Fontsize',labelSzAx);
ylim([0,pi*1.7]);
set(gca,'YTick',[0,pi],'YTickLabel',[0,180]);
xlabel('normalized distance');
% ylabel('\theta (deg)');
leg = legend('A-A''','B-B''','C-C''','D-D''');
set(leg,'Fontsize',labelLeg);
chld = get(gca,'Children');
set(chld,'Linewidth',2);
tx = text(0.18,pi*1.33,{...
    sprintf('\\sigma_{A-A''} = %0.0f^{\\circ}     \\sigma_{B-B''} = %0.0f^{\\circ}',std3a,std3b);...
    sprintf('\\sigma_{C-C''} = %0.0f^{\\circ}   \\sigma_{D-D''} = %0.0f^{\\circ}',std3c,std3d)});
    set(tx,'Fontsize',labelDeviations);
tx = text(0.05,pi*1.5,'f'); set(tx,'Fontsize',labelPlot);

pause(0.1);
subplot(3,2,1); set(gca,'Position',[xPos1 yPos1 wPos1 hPos]);
subplot(3,2,2); set(gca,'Position',[xPos2 yPos1 wPos2 hPos]);
subplot(3,2,3); set(gca,'Position',[xPos1 yPos2 wPos1 hPos]);
subplot(3,2,4); set(gca,'Position',[xPos2 yPos2 wPos2 hPos]);
set(leg,'Position',[xleg yleg wleg hleg]);
subplot(3,2,5); set(gca,'Position',[xPos1 yPos3 wPos1 hPos]);
subplot(3,2,6); set(gca,'Position',[xPos2 yPos3 wPos2 hPos]);




%% Saving

if saveArg == 1
    print(gcf,'cross_sections','-djpeg','-r300');
end

